close all;
clear all;

%% Block 0: Define Code Parameters 
% The quantification program functions to 
%   1) identify morphological regions
%   2) determine average pixel intensity
%   3) calculate fraction of pixels with elevated intensity

p_cut = 100; %Upper percentile for Endogenous Intensity Cutoff
blank_row = 2; %blank rows in directory
blank_row_training = 2; %blank rows in directory

%% Block 1A: Training Cutoff a-Syn
% To determine elevated pixel area we require determining the maximum
% endogenous control signal (cutoff_a and cutoff_p) for each experimental 
% replicate (Block 1), with a quality control parameter for high 
% autofluorescence images. Endogenous cutoff is determined by averaging the
% maximum pixel intensity of all contol images. Block 1A = a-Syn and Block
% 1B = PFF signal.

% Import Images
imagefiles_train_a = dir('/Users/adambindas/Dropbox/Documents/Koppes Lab/Programs/agg a-Syn Quantification/Att 4 01-23-21/Training Images (a-Syn)');

% Define Folder Size
nfiles_at = length(imagefiles_train_a)-blank_row_training;
training_cut_a = zeros(nfiles_at,1);

% Calculate Maximum Intensity for each Image
for ii=1:nfiles_at
   a_name = imagefiles_train_a(ii+blank_row_training).name; 
   a = imread(a_name);
   a_gray = rgb2gray(a);
   lower_thresh_a = prctile(a_gray,p_cut,'all');
   training_cut_a(ii)= lower_thresh_a;
   
end 

% Average the values for a single cutoff value.
cutoff_a = mean(training_cut_a)

%% Block 1B: Training Cutoff PFF

% Import Images
imagefiles_train_p = dir('/Users/adambindas/Dropbox/Documents/Koppes Lab/Programs/agg a-Syn Quantification/Att 4 01-23-21/Training Images (PFF)');

% Define Folder Size
nfiles_pt = length(imagefiles_train_p)-blank_row_training;
training_cut_p = zeros(nfiles_pt,1);

% Calculate Maximum Intensity for each Image
for ii=1:nfiles_pt
   p_name = imagefiles_train_p(ii+blank_row_training).name; 
   p = imread(p_name);
   p_gray = rgb2gray(p);
   lower_thresh_p = prctile(p_gray,p_cut,'all');
   training_cut_p(ii) = lower_thresh_p;
   
end 

% Average the values for a single cutoff value.
cutoff_p = mean(training_cut_p)

%% Block 1C: Importing Images and Preparing Variables
% This section of code imports remaining images and creates variables to be
% later filled through the remaining code.

% Import Image Files
imagefiles_asyn = dir('/Users/adambindas/Dropbox/Documents/Koppes Lab/Programs/agg a-Syn Quantification/Att 4 01-23-21/a-Syn');
imagefiles_pff = dir('/Users/adambindas/Dropbox/Documents/Koppes Lab/Programs/agg a-Syn Quantification/Att 4 01-23-21/PFF');
imagefiles_tubb = dir('/Users/adambindas/Dropbox/Documents/Koppes Lab/Programs/agg a-Syn Quantification/Att 4 01-23-21/B3T');
imagefiles_dapi = dir('/Users/adambindas/Dropbox/Documents/Koppes Lab/Programs/agg a-Syn Quantification/Att 4 01-23-21/DAPI');

% Define for loop length
nfiles = length(imagefiles_asyn)-blank_row;

% Create Variables for Quality Control
p_cut_low = 1; % Minimum pixel intensity
low_int = zeros(nfiles,1); % Fill variable for quality control

% Create remaining variables
asyn_avg_int = zeros(nfiles,1);
asyn_line_avg_int = zeros(nfiles,1);
asyn_soma_avg_int = zeros(nfiles,1);
asyn_tot_avg_int = zeros(nfiles,1);

asyn_p_area = zeros(nfiles,1);
asyn_p_area_line = zeros(nfiles,1);
asyn_p_area_soma = zeros(nfiles,1);
asyn_p_area_tot = zeros(nfiles,1);

pff_avg_int = zeros(nfiles,1);
pff_line_avg_int = zeros(nfiles,1);
pff_soma_avg_int = zeros(nfiles,1);
pff_tot_avg_int = zeros(nfiles,1);

pff_p_area = zeros(nfiles,1);
pff_p_area_line = zeros(nfiles,1);
pff_p_area_soma = zeros(nfiles,1);
pff_p_area_tot = zeros(nfiles,1);

cutoff_asyn = zeros(nfiles,1);
cutoff_asyn(:) = cutoff_a;
cutoff_pff  = zeros(nfiles,1);
cutoff_pff(:) = cutoff_p;

%% Block 2+: For Loop

for i = 1:nfiles % for loop to quantify each image

% Block 2a: Cutoff Intensity by Maximum Endogenous Signal
asyn_name = imagefiles_asyn(i+blank_row).name; % Extract file name
asyn = imread(asyn_name); % Open image

asyn_2 = rgb2gray(asyn); % Convert to grayscale
asyn_3 = asyn_2 - cutoff_a; % Subtract endogenous signal
asyn_4 = max(asyn_3,0); % Zero out all negative or non-elevate signal

a_im_low_int = prctile(asyn_2,p_cut_low,'all'); % Identify minimum intensity for quality control
low_int(ii) = a_im_low_int; % Fill minimum intensity

% Block 2b: PFF Cutoff + Comp Figure
pff_name = imagefiles_pff(i+blank_row).name; % Extract file name
pff = imread(pff_name); % Open image

pff_2= rgb2gray(pff); % Convert to grayscale
pff_3 = pff_2 - cutoff_p; % Subtract endogenous signal
pff_4 = max(pff_3,0); % Zero out all negative or non-elevate signal

% Block 3a: Mask Creation a-Syn
asyn_bw = im2bw (asyn_4,0.01); % Convert Elevated Pixels to a Black and White Image

% Block 3b: Mask Creation PFF
pff_bw= im2bw (pff_4,0.01); % Convert Elevated Pixels to a Black and White Image

% Block 3c: Quantification Prior to B3T Masks
%a-Syn Quantification
asyn_avg_int(i) = mean2(asyn_2); % Calculate average intensity of a-Syn for whole image
asyn_p_area(i) = sum(asyn_bw(:))/4194304; % Calculate fraction a-Syn elevated area for whole image

%PFF Quantification
pff_avg_int(i) = mean2(pff_2); % Calculate average intensity of a-Syn for whole image
pff_p_area(i) = sum(pff_bw(:))/4194304; % Calculate fraction a-Syn elevated area for whole image

% Block 4A: B3T Mask Neurites
tubb_name = imagefiles_tubb(i+blank_row).name; % Extract B3T filename
tubb = imread(tubb_name); % Import file name
tubb_2 = rgb2gray(tubb); % Convert to grayscale

background_tubb = imopen(tubb_2,strel('disk', 35)); % Identify background signal
tubb_3 = tubb_2 - background_tubb; % Subtract background signal
tubb_4 = max(tubb_3,0); % Zero negative values
tubb_5 = imadjust(tubb_4); % Saturate values
tubb_bw = im2bw (tubb_5,0.15); % Convert to black and white
tubb_bw_line = bwareaopen(tubb_bw,1000); % Remove small pixel count regions

% Block 4B: B3T Mask Somas
background_s_tubb = imopen(tubb_2,strel('disk', 500)); % Identify background signal
tubb_3s = tubb_2 - background_s_tubb; % Subtract background signal
tubb_4s = max(tubb_3s,0); % Zero negative values
tubb_5s = imadjust(tubb_4s); % Saturate values
tubb_bws = im2bw (tubb_5s,0.80); % Convert to black and white
tubb_bws2 = bwareaopen(tubb_bws,1500); % Remove small pixel count regions

%Dapi file
dap_name = imagefiles_dapi(i+blank_row).name; % Extract DAPI filename
dap = imread(dap_name); % Import file name
dap_2 = rgb2gray(dap); % Convert to grayscale

background_dap = imopen(dap_2,strel('disk', 500)); % Identify background signal
dap_3 = dap_2 - background_dap; % Subtract background signal
dap_4 = max(dap_3,0); % Zero negative values
dap_5 = imadjust(dap_4); % Saturate values
dap_bw = im2bw (dap_5,0.80); % Convert to black and white
dap_bw2 = bwareaopen(dap_bw,1500); % Remove small pixel count regions

tubb_bwt = max(tubb_bws2,dap_bw2); %Combine black and white DAPI and B3T masks
tubb_bw_soma = imfill(tubb_bwt,'holes'); % Fill in empty regions

% Block 4C: Combine B3T Mask and a-Syn/PFF Masks
tubb_bw_tot = max(tubb_bw_line,tubb_bw_soma); % Create a whole neuron mask

% Mask for Area Calculation
% Create a new b-w image of a-Syn elevated area within morphological regions
asyn_area_line = immultiply(asyn_bw,tubb_bw_line);
asyn_area_soma = immultiply(asyn_bw,tubb_bw_soma);
asyn_area_tot = immultiply(asyn_bw,tubb_bw_tot);

% Identify sum of a-Syn elevated pixels
asyn_mask_area_line = sum(asyn_area_line(:));
asyn_mask_area_soma = sum(asyn_area_soma(:));
asyn_mask_area_tot = sum(asyn_area_tot(:));

% Create a new b-w image of PFF elevated area within morphological regions
pff_area_line = immultiply(pff_bw,tubb_bw_line);
pff_area_soma = immultiply(pff_bw,tubb_bw_soma);
pff_area_tot = immultiply(pff_bw,tubb_bw_tot);

% Identify sum of PFF elevated pixels
pff_mask_area_line = sum(pff_area_line(:));
pff_mask_area_soma = sum(pff_area_soma(:));
pff_mask_area_tot = sum(pff_area_tot(:));

% Reduced Image for ROI Avg Int (restrict intensity vector to morphological
% region
asyn_roi_line = immultiply(imadd(asyn_2,1),tubb_bw_line);
asyn_roi_soma = immultiply(imadd(asyn_2,1),tubb_bw_soma);
asyn_roi_tot = immultiply(imadd(asyn_2,1),tubb_bw_tot);

pff_roi_line = immultiply(imadd(pff_2,1),tubb_bw_line);
pff_roi_soma = immultiply(imadd(pff_2,1),tubb_bw_soma);
pff_roi_tot = immultiply(imadd(pff_2,1),tubb_bw_tot);

% Block 5: Quantify Area
% a-Syn Area Quantification 
% Fraction of elevated pixels within morphological region
asyn_p_area_line(i) = sum(asyn_area_line(:))/sum(tubb_bw_line(:));
asyn_p_area_soma(i) = sum(asyn_area_soma(:))/sum(tubb_bw_soma(:));
asyn_p_area_tot(i) = sum(asyn_area_tot(:))/sum(tubb_bw_tot(:));

% Average intensity within morphological region
asyn_line_avg_int(i) = mean2(nonzeros(asyn_roi_line))-1;
asyn_soma_avg_int(i) = mean2(nonzeros(asyn_roi_soma))-1;
asyn_tot_avg_int(i) = mean2(nonzeros(asyn_roi_tot))-1;

%PFF Quantification
% Fraction of elevated pixels within morphological region
pff_p_area_line(i) = sum(pff_area_line(:))/sum(tubb_bw_line(:));
pff_p_area_soma(i) = sum(pff_area_soma(:))/sum(tubb_bw_soma(:));
pff_p_area_tot(i) = sum(pff_area_tot(:))/sum(tubb_bw_tot(:));

% Average intensity within morphological region
pff_line_avg_int(i) = mean2(nonzeros(pff_roi_line))-1;
pff_soma_avg_int(i) = mean2(nonzeros(pff_roi_soma))-1;
pff_tot_avg_int(i) = mean2(nonzeros(pff_roi_tot))-1;

% Report progress of for loop
i
nfiles

end

%% Block 6 Make into a table
% Extract row names
rownames = transpose({imagefiles_asyn(blank_row+1:blank_row+nfiles).name});

% baseline intensity quality issue identification
% low_int_median = median(low_int);
% low_int_qual = low_int > low_int_median+2;
% info_qual = horzcat(low_int,low_int_qual);

column_names = {
    'rownames',
    'cutoff_asyn',
    'cutoff_pff',
    'asyn_avg_int', 
    'asyn_roi_line_avg_int', 
    'asyn_roi_soma_avg_int', 
    'asyn_roi_tot_avg_int',
    'asyn_p_area',
    'asyn_p_area_line',
    'asyn_p_area_soma',
    'asyn_p_area_tot',
     'pff_avg_int', 
     'pff_roi_line_avg_int', 
     'pff_roi_soma_avg_int', 
     'pff_roi_tot_avg_int',
     'pff_p_area',
     'pff_p_area_line',
     'pff_p_area_soma',
     'pff_p_area_tot'
    };

% att_csv = horzcat(rownames,asyn_avg_int,asyn_line_avg_int,asyn_soma_avg_int,asyn_tot_avg_int,asyn_p_area,asyn_p_area_line,asyn_p_area_soma,asyn_p_area_tot,pff_avg_int,pff_line_avg_int,pff_soma_avg_int,pff_tot_avg_int,pff_p_area,pff_p_area_line,pff_p_area_soma,pff_p_area_tot);
att_noname_csv = table(rownames,cutoff_asyn,cutoff_pff,asyn_avg_int,asyn_line_avg_int,asyn_soma_avg_int,asyn_tot_avg_int,asyn_p_area,asyn_p_area_line,asyn_p_area_soma,asyn_p_area_tot,pff_avg_int,pff_line_avg_int,pff_soma_avg_int,pff_tot_avg_int,pff_p_area,pff_p_area_line,pff_p_area_soma,pff_p_area_tot);
imaging_data_csv = table(rownames,cutoff_asyn,cutoff_pff,asyn_avg_int,asyn_line_avg_int,asyn_soma_avg_int,asyn_tot_avg_int,asyn_p_area,asyn_p_area_line,asyn_p_area_soma,asyn_p_area_tot,pff_avg_int,pff_line_avg_int,pff_soma_avg_int,pff_tot_avg_int,pff_p_area,pff_p_area_line,pff_p_area_soma,pff_p_area_tot,'VariableNames',{'rownames','cutoff_asyn','cutoff_pff','asyn_avg_int', 'asyn_roi_line_avg_int', 'asyn_roi_soma_avg_int', 'asyn_roi_tot_avg_int','asyn_p_area','asyn_p_area_line','asyn_p_area_soma','asyn_p_area_tot','pff_avg_int', 'pff_roi_line_avg_int', 'pff_roi_soma_avg_int', 'pff_roi_tot_avg_int','pff_p_area','pff_p_area_line','pff_p_area_soma','pff_p_area_tot'});

writetable(imaging_data_csv)
